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13.  The  major  goal  to  this  research  was  to  develop  a  strategy  for  establishing  a 
microscopic  atomic-level  understanding  of  the  fundamental  surface  processes 
ultimately  responsible  for  friction,  adhesion  at  surfaces,  and  abrasion.  The 
approach  was:  (i)  to  use  quantum  chemical  studies  to  establish  the  dominant 
surface  species  for  clusters  of  atoms  modeling  various  ceramics  and  to 
elucidate  tna  thermochemistry  and  detailed  mechanistic  steps  involved  in 
surface  reactions  of  such  systems;  (ii)  to  develop  theoretical  force  fields 
based  on  the  energy  surfaces  from  clusters  in  i  that  allow  predictions  of  the 
energies  and  geometries  for  infinite  surfaces  and  interfaces;  (iii)  to  use  the 
force  fields  from  ii  to  predict  the  barriers  and  kinetics  for  various 

diffusion  and  reaction  processes  relevant  for  catalysis,  corrosion,  and 
materials  synthesis  processes;  (iv)  to  develop  procedures  for  molecular 
dynamics  and  Monte  Carlo  simulations  of  various  chemical  processes  in  these 
systems,  and,  (v)  to  interface  the  results  of  these  simvlati  or'1'  onto 
appropriate  graphics  systems,  allowing  the  designer  to  interactively  follow  a 
three-dimensional  image  of  the  evolving  system. 

Through  simultaneous  examination  of  all  properties  of  the  surface  < chemical  . 
tribological,  physical,  and  mechanical)  with  the  same  theory,  we  have  the 

possibility  of  strong  tests  on  different  aspects  of  the  theory.  In  addition, 
such  a  central  theory  would  provide  new  connections  between  these  properties 
that  would  serve  to  connect  what  are  now  orthogonal  experiments.  These 
developments  should  provide  a  level  of  understanding  useful  in  designing  new 
materials  and  such  computer-aided  materials  simulations  (CAMS)  should  provide 
tools  allowing  many  of  the  design  concepts  to  be  tested  on  the  computer  in 

advance  to  attempting  difficult  syntheses  and  characterizations  in  the 

laboratory . 
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SUMMARY 

The  goal  of  this  research  was  to  develop  a  microscopic  atomic-level  understanding  of 
the  fundamental  surface  processes  responsible  for  friction,  adhesion,  and  abrasion.  The 
emphasis  is  on  materials  (particularly  ceramics)  having  potential  applications  in  high  per¬ 
formance  systems  and  also  on  systems  for  which  it  is  possible  to  design  precise  experimental 
tests  of  the  predictions. 

There  have  been  five  major  advances: 

a.  Pseudospectral-Generalized  Valence  Bond  (PS-GVB)  for  accurate  quantum 
chemical  calculations  of  large  molecules.  This  is  a  major  breakthrough  for  calculating 
electronic  states  and  properties  of  clusters  and  crystals.  With  previous  methods,  in¬ 
clusion  of  electron  correlation  effects  (as  in  GVB)  has  generally  been  limited  to  about 
10  atoms  or  so  (200  basis  functions)  where  costs  scales  as  N 4  or  Ns  with  size  (N) 
of  this  system.  PS-GVB  avoids  calculating  the  Ar4  electron  repulsion  integrals  which 
causes  this  problem.  Instead  the  Coulomb  and  exchange  operators  are  calculated  di¬ 
rectly  (over  a  numerical  grid)  and  then  combined  to  form  the  Foci:  operators  from 
which  the  GVB  orbitals  are  obtained.  This  decreases  the  computational  orobiem  to 
one  that  scales  as  V'3  (dropping  to  A'3  with  cutoffs)  and  allows  systems  of  2000  basis 
functions  100-200  atoms.'  to  be  commuted  in  the  same  time  as  nrcviousiv  for  200 
basis  functions  ,10-20  atoms  i.  This  approach  also  allows  use  of  periodic  bc-undarv 
conditions  (PBC)  for  direct  calculation  of  infinite  crystals  and  surfaces.  (PBC  is  still 
under  development.) 

b.  The  Cell  Alultipole  Alethod  (CMM)  for  calculating  forces,  structures,  and  proper¬ 
ties  of  million-atom  systems.  This  is  a  major  breakthrough  for  atomic-level  molecular 
dynamics  simulations  of  condensed  systems.  With  previous  methods,  simulations  on 
thousands  of  atoms  generally  requires  CRAY  level  computers  for  studying  dynamical 
properties;  however,  for  modeling  tribiological  systems  it  may  be  necessary  to  treat 
100,000  or  even  a  million  atoms.  (A  diamond  cube  of  200  A  —  20  nm  =  0.02  q  per  side 
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has  one  million  atoms.)  Although  there  are  N2  Coulomb  interactions  for  an  N  body 
systems  the  computational  cost  for  the  Cell  Multipole  Method  is  rigorously  linear  in 
N.  This  is  done  by  evaluating  multipole  fields  for  each  region  (cell).  We  have  carried 
out  million  atom  simulations  for  finite  polymers  and  for  million  atoms  per  unit  cell 
simulations  on  pseudo  periodic  amorphous  systems.  In  both  cases  CMM  is  103  faster 
than  the  explicit  evaluation  of  all  interactions. 

c.  The  Charge  Equilibration  (QEq)  approach  for  predicting  charge  distributions  in 
ceramic  and  polymer  systems.  This  approach  solves  a  major  problem  in  simulation*  cf 
ceramics  and  polymer  systems  (where  electrostatic  effects  play  a  dominant  role)  since 
th^re  had  been  no  method  for  actually  predicting  the  charge  distributions  essential  to 
simulations.  The  new  method  has  as  its  input  only  atomic  information  (available  for 
all  atoms)  and  predicts  the  charge  redistribution  as  a  function  of  geometry. 

d.  The  Hessian-Biased  Force  Field  (HBFF)  method  of  combining  data  (the  Hessian 
or  second-derivative  matrix)  from  quantum  chemical  calculations  with  observed  vi¬ 
brational  frequencies  to  determine  accurate  valence  force  fields  (describing  how  the 
energy  and  forces  depend  on  bonds,  angles,  etc.).  This  is  critical  for  describing  the 
dynamics  of  motion  as  two  surfaces  interact  with  each  other  or  with  lubricant,  and 
it  is  also  essential  in  describing  the  vibrational  modes  (phonon  bands)  that  couple  to 
surface  states  in  dissipative  (frictional)  processes.  These  data  are  also  essential  for 
calculating  quantities  related  to  the  hardness  of  materials. 

e.  The  Interstital  Electron  Model  (IEM)  for  Simulations  of  Metallic  Systems. 
We  believe  that  this  is  a  major  breakthrough  that  will  allow  meaningful  simulations 
of  tribological  properties  of  metals.  The  difficulty  in  simulating  metals  is  that  the 
electrons  respond  dramatically  to  motion  of  the  ions  and  hence  the  forces  cannot  be 
described  solely  in  terms  of  ion  positions.  Our  new  approach  (based  on  concepts  arising 
from  GVB  quantum  mechanical  calculations)  treats  electrons  as  classical  particles  so 
that  the  energy  is  described  in  terms  of  electron-electron,  electron-ion,  and  ion-ion 
interactions.  We  have  tested  this  approach  by  calculating  elastic  constants  and  phonon 
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dispersion  curves,  with  excellent  results. 

In  addition  we  have  developed  force  fields  suitable  for  simulations  on  graphite,  dia¬ 
mond,  fullerenes  (C6o,  C70  buckyballs),  Sij  N4,  SiC,  Si,  Ge,  Sn,  cubic-BN  GaAs,  GaP, 
GaSb,  AlAs,  A1P,  AlSb,  InAs,  InP,  InSb,  SiC>2,  AI2  Oj,  NiO,  YBa2Cu30a:,  Al,  Fe,  Ni,  Cu, 
Pd,  Ag,  Pt,  Au,  and  other  materials  of  tribiological  interests. 

We  believe  that  these  tools  and  force  fields  provide  the  basis  for  microscopic  theoretical 
modeling  of  the  Chemical  and  Tribiological  Properties  of  Ceramic  Surfaces  and  Interfaces. 

I.  Introduction 

The  major  goal  to  this  research  was  to  develop  a  strategy  for  establishing  a  microscopic 
atomic-level  understanding  of  the  fundamental  surface  processes  ultimately  responsible  for 
friction,  adhesion  at  surfaces,  and  abrasion. 

Although  there  have  been  rapid  advances  in  the  experimental  techniques  for  examin¬ 
ing  composition  at  surfaces  and  interfaces  relevant  for  catalysis,  tribology,  corrosion,  and 
materials  synthesis,  there  has  been  very  little  in  the  way  of  a  microscopic  theoretical 
model  suitable  for  understanding  the  chemical,  physical,  and  mechanical  properties  in 
terms  of  atomic-level  structure  and  bonding  concepts.  Our  project  was  a  broad-based 
theoretical  program  aimed  at  providing  the  fundamentals  for  constructing  such  an  atomic- 
level  microscopic  model  of  the  chemical,  tribological ,  physical,  and  mechanical  properties  of 
ceramic  surfaces  and  interfaces . 

The  approach  was 

:.  to  use  quantum  chemical  studies  to  establish  the  dominant  surface  species  for  clusters 
of  atoms  modeling  various  ceramics  and  to  elucidate  the  thermochemistry  and  detailed 
mechanistic  steps  involved  in  surface  reactions  of  such  systems; 
ii.  to  develop  theoretical  force  fields  based  on  the  energy  surfaces  from  clusters  in  :  that 
allow  predictions  of  the  energies  and  geometries  for  infinite  surfaces  and  interfaces; 

1'::.  to  use  the  force  fields  from  ii  to  predict  the  barriers  and  kinetics  for  various  diffu¬ 
sion  and  reaction  processes  relevant  for  catalysis,  corrosion,  and  materials  synthesis 
processes; 
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it),  to  develop  procedures  for  molecular  dynamics  and  Monte  Carlo  simulations  of  various 
chemical  processes  in  these  systems;  and, 

v.  to  interface  the  results  of  these  simulations  onto  appropriate  graphics  systems,  allowing 
the  designer  to  interactively  follow  a  three-dimensional  image  of  the  evolving  system. 
Through  simultaneous  examination  of  all  properties  of  the  surface  (chemical,  tribological, 
physical,  and  mechanical)  with  the  same  theory,  we  have  the  possibility  of  strong  tests 
on  different  aspects  of  the  theory.  In  addition,  such  a  central  theory  would  provide  new 
connections  between  these  properties  that  would  serve  to  connect  what  are  now  orthogonal 
experiments.  These  developments  should  provide  a  level  of  understanding  useful  in  design¬ 
ing  new  materials  and  such  computer-aided  materials  simulations  (CAMS)  should  provide 
tools  allowing  many  of  the  design  concepts  to  be  tested  on  the  computer  in  advance  to 
attempting  difficult  syntheses  and  characterizations  in  the  laboratory. 

There  have  been  major  advances  in  each  of  these  areas,  as  is  summarized  in  the  next 
few  sections. 

II.  New  Developments  in  Quantum  Chemistry 

Over  the  last  ten  years,  the  Goddard  group  has  carried  out  a  number  of  quantum 
chemical  studies  on  clusters  designed  to  mimic  the  chemically  active  sites  of  metal,  semi¬ 
conductor,  and  ceramic  surfaces.  This  research  has  led  to  the  GYB  suite  of  program.' 
which  provides  powerful  methods  of  including  electron-correlation  (many-body)  effects  for 
accurate  predictions  of  reaction  intermediates  in  transition  metal  complexes.  These  studies 
have  led  to  a  number  of  insights  on  surface  chemistry  and  have  led  to  new  or  modified 
mechanisms  for  key  steps  of  several  important  catalytic  reactions. 

Despite  the  progress,  there  remained  severe  limitations  for  use  of  these  techniques  and 
tribological  systems.  We  are  convinced  that  ab  initio  methods  including  a  high  level  of  elec¬ 
tron  correlation  are  essential  to  obtain  accurate  information  about  the  surface  chemistry 
involved  in  wear  and  friction  (since  there  is  little  experimental  data  to  use  as  guidelines  in 
adjusting  parameters).  However,  with  previous  techniques  inclusive  of  electron  correlation 
has  restricted  use  to  using  small  clusters  as  models  of  catalytic  systems,  making  it  difficult 
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to  examine  the  role  of  promoters,  poisons,  modifiers,  solvents,  and  surface  structures  on 
catalytic  pathways  and  rates. 

A  fundamental  advance  here  has  been  the  development  of  PS-GVB  (pseudospectral- 
generalized  valence  bond)  where  the  Coulomb  and  exchange  potentials  are  calculated  di¬ 
rectly  (over  a  numerical  grid)  without  evaluating  or  storing  the  two-electron  integrals. 
Since  in  PS-GVB  the  N 4  set  of  two-electron  integrals  are  never  constructed,  we  eliminate 
the  source  of  bottlenecks  in  previous  ab  initio  quantum  chemistry  methods,  storing  and 
manipulating  N4  of  these  integrals.  Thus,  for  200  basis  functions  the  usual  approach  in¬ 
volved  AM/8  =  2  x  10®  integrals;  for  2000  basis  functions  there  are  2  x  1 0 1 2 !  With  PS,  the 
Coulomb  and  exchange  integrals  are  evaluated  over  a  numerical  grid  and  used  as  needed. 
For  HF  wavefunctions  with  200  basis  functions,  the  PS-HF  program  on  a  CRAY-YMP  is 
10  times  faster  than  GAUSSIAN  88  (20  times  faster  than  GAUSSIAN  86),  and  for  2000 
basis  functions  it  should  be  100  times  faster. 

There  are  four  points  to  the  overall  PS-GVB  development: 
i.  PS-GVB-PP  allows  calculation  of  structures  and  charges  for  large  clusters  (20-200 
atoms)  useful  in  establishing  parameters  for  modeling  realistic  materials  problems. 
GVB-PBC  will  provide  correlated  wavefunctions  j’PP)  for  infinite  periodic  crystals 
(up  io  ~50  adorns  per  ceil;,  allowing  direct  prediction  of  the  electronic  and  structural 
properties  for  crystals  and  surfaces  (without  corrections  due  to  cluster  effects!. 

:::.  ?S-G\  B-R.Ci  will  provide  accurate  descriptions  of  saddle  noir.ts  for  reactions  i  up  to 
—  30  atoms;,  arrowing  redtstic  simulation  of  ca*alysts  fmd  rear*; on 


G\B-EN2  (second  order  Epstein-Nesbet  perturbation  theory  for  GVB-R.CI  wavefunc¬ 
tions)  should  allow  very  accurate  description  of  bond  energies  (~0.5  kcal/mol)  and 
van  der  V  aals  interactions  (to  obtain  nonbonded  potentials  for  molecular  dynamics 
studies). 

Ve  consider  PS-GV  B  to  be  a  geniune  breakthrough,  providing  the  means  for  practical 
calculations  from  first  principles  the  chemical  and  mechanical  properties  of  surfaces  and 
interfaces. 


These  PS-GVB  methods  are  under  active  development,  and  all  are  expected  to  be 


fully  implemented  by  summer  1992. 

A.  Pseudospectral  Methods 

In  PS-GVB  the  usual  N4/S  set  of  two-clectron  integrals  (/j.u\arj)  is  never  evaluated 
(where  N  is  the  number  of  basis  functions).  It  is  storing  and  manipulation  of  these  inte¬ 
grals  that  causes  the  current  bottlenecks  in  quantum  chemical  calculations  (for  200  basis 
functions  there  are  N4  / 8  =  2  x  10s  integrals,  for  2000  basis  functions  there  are  2  *  1012). 
Instead,  with  PS  the  Coulomb  and  exchange  operators  are  evaluated  directly  for  a  numer¬ 
ical  grid.  This  involves  evaluating  the  TV2/ 2  quantities 
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directly  for  a  set  of  grid  points  {rfl}.  With  Gaussian  type  basis  functions,  the  potentials 
Acrr]{Tg)  are  evaluated  analytically  for  each  grid  point  (a  sufficiently  fast  calculation  that 
these  potentials  can  be  used  and  then  discarded). 

Given  the  potentials  (l),  the  fundamental  Coulomb  and  exchange  matrices  could  be 

calculated  numerically  using 

J^u  =<  m!  JjW  >=  xMxM  ^  c<rjcijA^v [9) 

9  c 9  ,ns 

=<  ^\Kj\v  >=  X]  X^g)  L  '  H  CejC^jArvis) 


no weve:  2)  ieacs  to  tame  numeric 
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JjX*  and  KjXv  (3) 

have  components  outside  the  basis  space  ^Xu}-  ^  hmte  grids  this  leads  to  large  an  using 
errors,1  in  (2)  which  manifests  itself  as  noise  in  the  wavefunction.2  The  solution  to  this 
difficulty  is  to  expand  (3)  in  terms  of  a  larger  deahasing  basis  {x*;c  —  1  ...M’} 

J:xM  = 
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used  only  for  evaluating  (5) 
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where  the  5*^.  =<  Xu Ix^-  >  the  overlap  between  functions  of  the  normal  and  augmented 
(dealiasing)  basis  sets.  This  set  of  dealiasing  functions  filters  out  the  noise  due  to  incom¬ 
pleteness  of  {x^}-  The  coefficients  and  in  (4)  are  obtained  by  a  least  squares 

fit  over  the  grid  of  points.3 

Critical  to  the  efficiency  of  the  PS  method  is  the  multigrid  strategy  (involving  medium, 
fine,  and  ultrafine  grids)  that  is  used  in  different  parts  of  the  calculation.  For  HF  wave- 
functions  these  procedures  have  been  well  optimized  by  Friesner  e£  a/.3  and  for  molecules 
with  200  basis  functions  the  current  PS-HF  program  on  a  CRAY-YMP  is  10  times  faster 
than  GAUSSIAN  88  (20  times  faster  than  GAUSSIAN  86).  For  systems  above  this  size 
computation  times  for  Gaussian  88  scales  as  A'3  whereas  PS-HF  scales  as  A’2.  Thus  for 
2000  basis  functions  (~  100  atoms)  PS-EF  should  be  100  times  faster.  Note  that  2000 
basis  functions  makes  sense  for  PS  since  storage  and  IO  needs  go  as  A  "  (A'*  with  cut  off 
rather  than  A  "  (A  3  with  cutoffs). 

The  Goddard  croup  in  collaboration  with  the  Friesner  croup  has  shown  that  the  PS 

.vcfunctions.  where  the  "ouiomD  and  cncnar.rv 
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are  usee  direct,  v 
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G Y3-PP  program  is  operational  and  has  been  used  for  tests  on  singlet-triplet  excitation 
energies,  rotational  barriers  about  double  bonds,  and  bond  dissociation. ’  In  addition,  it 
has  been  tested  for  use  with  effective  potentials.'1'5  In  all  these  cases  we  ffnd  better  than  0.1 
kcal/’mole  absolute  accuracy  and  0.01  kcal/mole  relative  accuracy,  as  compared  to  standard 
approaches. 

We  are  in  the  process  of  completely  redeveloping  both  the  GYB  formulation  and 
the  PS  techniques,  grids  and  cutoffs  for  correlated  pairs.  Although  a  challenging  goal, 
we  believe  that  calculations  with  2000  basis  functions  and  100  GYB  pairs  will  become 


9 


possible.  Tliis  would  allow  accurate  studies  of  ceramic  clusters  with  up  to  200  atoms  and 
transition  metal  systems  with  up  to  100  atoms,  and  would  allow  accurate  calculations  for 
complexes  providing  good  models  of  important  surface  processes. 

An  important  virtue  of  the  PS  approach  is  that  the  exchange  terms  needed  for  G\  B- 
RCI  wavefunctions  can  be  calculated  directly.  The  GYB-RCI  wavefunction6’7  (which  in¬ 
cludes  all  spin  coupling  and  interpair  correction  terms)  leads  to  quite  accurate  potential 
curves  for  reactions  i/ the  orbitals  are  optimized  self  consistently.  (In  this  case,  the  wave- 
function  contains  essentially  all  important  parts  of  a  full  GVB-CI  or  CAS-SGF  wavefunc¬ 
tion.)  However  in  the  past  such  self  consistent  GVB-RCI  calculations  required  use  of  the 
GYB3  program  and  hence  an  A75  integral  transformation  every  iteration.  Fortunately 
with  PS  we  can  evaluate  the  additional  exchange  potentials  required  in  GYB-RCI  directly 
(on  the  fly  without  calculating  two-electron  integrals).  Thus  with  PS  we  can  generalize 
the  GYB-PP  program  to  calculate  the  GVB-RCI  wavefunction  self  consistently  without 
integral  transformations.  This  requires  a  modification  in  how  spin  couplings  are  treated, 
with  a  slight  restriction  on  the  generality  of  the  spinfunctions.  Thus  for  such  very  ac¬ 
curate  wavefunctions,  the  costs  should  drop  from  A5  to  A73  and  with  cutoffs  to  A"3,  for 


lame  svstems.  Timing  estimates  are  more  difficult  here,  but  we  believe  that  sel 
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-GYB-R.CI  for  up  to  1000  basis  functions  and  50  GY3  pairs  might  be  possible 


are  currently  d; 


velcoir.c  the  theory  for  PS-GVB-RCI  and  e::u>c 
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For  tribological  orobiems.  it  is  essential  to  use  nertodic  boundary  conditions  (PBC)  in 


ing  the  wavefunctions  (obviating  the  need  for  cluster  approximations).  For  closed 


sne.i  rrr  wav 


efunctior.s  this  is  straightforward,  but  the  local  exchange  terms  converge  slowly 


with  distance,  leading  to  long  expansions.9  By  taking  suitable  Fourier  transformations,  it  is 
possible  to  do  the  slowly  convergent  integrals  in  reciprocal  space10;  however  this  requires  a 
complicated  integral  transformation.  With  PS  this  PBC  process  is  considerably  simplified 
because  (1)  can  be  used  to  calculate  the  exchange  terms  as  they  are  needed.  Thus  using 
PS  it  is  quite  practical  to  calculate  HF  wavefunctions  using  periodic  boundary  conditions 
(KF-PBC).  tor  GYB  wavefunctions  use  of  PBC  is  more  complicated  because  the  individual 
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GVB  Fock  operators  do  not  combine  into  a  single  effective  multiparticle  operator.  However 
for  completely  filled  bands  it  is  straightforward  to  reformulate  GVB-PP  wavefunctions  for 
PBC.  This  requires  additional  localized  exchange  terms,  but  these  can  be  calculated  using 
the  PS  approach.  The  restriction  to  closed  shells  (filled  bands)  will  be  no  limitation  for 
many  cases  (semiconductors  and  ceramics),  although  it  precludes  study  of  metallic  systems 
(where  partially  filled  bands  are  necessary).  GVB-PBC  will  allow  us  to  treat  the  full  crystal 
for  many  interesting  catalysts.  We  are  currently  developing  algorithms  for  the  GVB-PBC 
program11  and  expect  to  have  a  working  version  by  early  1992.  Timing  estimates  are 
premature  but  our  goal  is  to  eventually  treat  systems  with  1000  basis  functions  and  50 
GVB  pairs  per  unit  cell  (50  transition  atoms  per  unit  cell). 

Such  calculations  of  three-dimensional  crystals  should  be  very  valuable  in  extracting 
charges  and  force  fields  that  could  be  used  for  molecular  mechanics  and  molecular  dynamics 
calculations  of  tribological  systems.  By  using  finite  thickness  slabs  we  will  use  PBC  for 
crystal  surfaces,  allowing  us  to  consider  structures  and  energetics  for  various  reaction 
intermediates.  (Since  PBC  is  used,  the  fundamental  repeating  unit  must  be  large  enough 
that  the  equivalent  catalytic  sites  do  not  interact.)  For  GVB-RCI  a  general  formulation 
of  PBC  also  seems  possible,  with  the  restriction  to  filled  bands,  but  we  have  not  worked  it 
out.  However,  including  the  RCI  spin  couplings  for  the  two  to  four  pairs  that  recouple  in 
the  transition  state  should  be  possible. 

With  PS  it  is  possible  to  start  with  GVB-RCI  wavefunctions  and  calculate  Epstein- 
Mesbet  perturbation  theory  though  second  order.  This  should  allow  bond  energies  accurate 
to  ~0.5  kcal/mol  and  should  allow  ab  initio  determination  of  parameters  for  van  der  Waals 
(dispersion)  interactions  important  for  treating  lubricants  interfacing  with  solid  surfaces. 
The  theory  has  been  partly  worked  out  and  the  initial  GVB-RCI-EN2  programs  should  be 
developed  by  the  Summer  of  1992. 

III.  New  Developments  in  Force  Fields 

Despite  the  progress  in  treating  larger  clusters  with  ab  initio  quantum  chemistry,  the 
calculations  are  far  too  slow  for  studying  the  dynamics  of  realistic  models  of  catalysts.  Thus 
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it  is  essential  to  replace  the  electrons  with  a  force  field  suitable  for  molecular  dynamcis 
and  Monte  Carlo  simulations.  Typically12  force  fields  are  described  in  terms  of  short-range 
valence  (bonded)  interactions  ( Evai )  expressed  in  terms  of  two-body  ( Ebond )>  three-body 
(■^cnj/c))  and  four-body  i^&toraion-)  ^inversion')  terms, 


E val  —  Ebond  -f-  Ean„[e  -f-  Etoraion  "h  E{ 


nvcrMxorx 


(6a) 


and  long-range  nonbonded  interactions  (-Enb)  composed  of  van  der  Waals  (Evcjw)  and 
electrostatic  (Eq)  terms 

Ent,  =  E„d  w  +  Eq  (6  b) 

A.  Electrostatics  (Charge  Equilibration) 

Given  a  set  of  charges  Q{  on  various  atoms,  the  electrostatic  interaction  is 


£«  =  £ 

3 


QiQj 

eRij 


Traditionally,  in  molecular  simulations  of  organic  molecules,  the  charges  are  estimated 
based  on  Hartree-Fock  (HF)  calculations  (with  the  charges  chosen  to  £t  the  calculated 
potential  energy  distribution).  There  are  several  serious  problems  with  this  approach: 

c.  Accurate  HF  calculations  (631G*X  basis)  are  possible  only  for  small  molecules,  un  to 
10-20  atoms  or  so  (with  PS-HF  it  should  become  possible  for  100  atoms  A 

b.  The  charges  should  be  allowed  to  change  as  atoms  move  about,  as  bonds  are  broken, 
etc.,  but  HF  calculations  every  time  step  is  not  practical. 

c.  In  order  to  recalculate  the  charges  during  the  dynamics,  the  charge  calculation  must 
be  extremely  rapid  and  cannot  involve  calculating  wavefunctions. 

d.  There  are  no  procedures  for  estimating  the  charges  for  crystals  of  ceramics,  metals, 
etc.  important  to  tribiology. 

M'hat  is  needed  is  a  general  approach  that  can  provide  charges  for  large  or  infinite  systems 
while  allowing  the  charges  to  polarize  and  change  during  the  dynamics.  We  believe  that 
we  have  solved  this  critical  problem  in  a  recent  series  of  papers.13 
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This  new  method  is  called  Charge  Equilibration  (QEq),  and  the  basic  idea  is  to 
consider  each  atom  as  a  separate  system  in  a  grand  canonical  ensemble.  The  energy  of 
each  atom  (A)  is  considered  to  depend  quadratically  upon  the  charge  Q  a 

Ea(Q  a)  =  Eq  +  tjaQ  a  +  \J  aaQ\  (8) 

where  tja  and  J a  are  related  to  electronegativity  and  atomic  Coulomb  repulsion  (electrical 
capacity),  respectively.  The  interaction  between  two  atoms  is  written  in  the  form 

Jab{R)QaQb  (9) 

where  the  interaction  term  Jab(R)  has  the  Coulomb  form  (I/R^b)  f°r  large  R,  but  in¬ 
cludes  shielding  so  that  it  goes  to  a  finite  limit  as  R  — *  0.  Requiring  that  the  charge 
distribution  be  in  equilibrium  (constant  chemical  potential)  leads  then  to  a  set  of  linear 
equations 

AijQj  =  Bi  (10) 

fox  the  equilibrium  charges  Q  at  a  given  geometry. 

The  QEq  method  is  formulated  so  that  there  are  no  free  parameters.  The  only  data 
are  the  experimental  ionization  potential  (IP)  and  electron  affinity  (EA)  (which  determine 
t] a  and  J aa)  and  the  atomic  radii  [which  determine  the  shielding  of  J ab{R )  for  small  R. i. 
The  resulting  charges  are  in  good  agreement  with  those  obtained  by  fitting  electrostatic 
potentials  from  good  quality  Hr  and  MP2  wavefunctions.  This  represents,  we  believe,  a 
significant  advance  in  simulations. 

The  QEq  approach  requires  quantitative  electronegativity  and  self  Coulomb  terms 
from  experimental  data.  We  found  a  way  to  generalize  ideas  of  Mulliken  and  Pauling  to 
obtain  an  unambiguous,  quantitative  measure  of  these  quantities  directly  from  experimen¬ 
tal  atomic  data.13  This  approach  is  useful  for  both  main  group  elements  and  transition 
metals.13 

The  Charge  Equilibration  method  (QEq)  uses  only  experimental  atomic  data  (ion¬ 
ization  potential,  electron  affinity,  and  bond  radius)  and  can  be  used  for  anjr  atoms  up 
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through  Lr  (element  103).  With  current  programs,  molecules  containing  600  atoms  are 
quite  practical  (there  is  no  hard  limit).  This  QEq  method  has  now  been  extended  to  in¬ 
finite  systems  (PBC)  where  crystals  with  200  atoms  per  unit  cell  are  quite  practical  (no 
hard  limit).  In  addition,  the  dynamic  charge  fluctuations  during  vibrations  have  recently 
been  included.  This  is  critical  for  describing  the  optical  vibrational  modes  of  nonmetals. 
In  addition  these  dynamical  changes  will  allow  infrared  and  Raman  intensities  to  be  pre¬ 
dicted.  These  methods  have  been  extended  to  allow  prediction  of  dielectric  constants  and 
piezoelectric  constants. 

In  addition  to  the  above  semiempirical  approach  for  charges,  we  have  been  developing 
a  first  principles  method  suitable  for  large  molecules  and  crystals  (up  to  200  atoms).  Using 
the  PS  approach  we  have  now  modified  the  PS-GVB-PP  program  to  calculate  the  charges 
based  on  electrostatic  potentials  of  correlated  wavefunctions.  [Up  to  2000  basis  functions 
200  atoms).] 


B.  Valence  Force  Fields 

The  most  common  approach  to  developing  valence  force  fields  for  molecular  simula¬ 
tions  has  been  to  start  with  the  experimental  geometries  and  vibrational  frequencies  for 
a  set  of  related  molecules  and  to  find  force  field  parameters  that  provide  the  best  fit  to 
these  data.  Unfortunately,  the  experimental  vibrational  data  do  not  provide  sufficient 
constraints  to  ensure  that  the  second  derivative  matrix  (Hessian)  is  well  fit.  To  illustrate 
the  problem,  consider  a  general  nine-atom  molecule.  A  simple  spectroscopic  quality  force 
field  would  involve  a  minimum  of  18  parameters  to  describe  the  nine  different  bond  stretch 
terms,  104  parameters  to  describe  the  18  different  bond  angle  terms,  and  at  least  24  pa¬ 
rameters  to  describe  the  24  different  torsion  terms.  Thus  146  parameters  are  needed  for 
the  force  field.  However,  the  21  normal  modes  provide  only  21  pieces  of  data  to  determine 


second  derivatives,  and  the  experimental  structure  provides  only  21  constraints  on  the  first 
derivatives  (forces).  These  42  experimental  constraints  are  far  short  of  determining  the 
146  FF  parameters.  Our  conclusion  is  that  to  generate  accurate  force  fields,  one  must  use 
additional  information  beyond  that  available  from  spectroscopic  studies. 
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An  alternative  is  to  use  the  power  of  modern  quantum  chemistry  to  calculate  the  full 
second  derivative  matrix  (Hessian)  at  any  geometry  for  simple  wavefunctions.  Unfortu¬ 
nately,  simple  HF  wavefunctions  lead  to  vibrational  frequencies  that  typically  disagree  with 
experiment  by  ~  10  —  20%.  One  can  obtain  greater  accuracy  with  correlated  wavefunctions 
(MP2,  GVB,  Cl);  however,  full  Hessians  from  such  wavefunctions  are  currently  practical 
only  for  tiny  systems.  [The  development  of  PS-MP2  and  PS-GVB-RCI  may  change  this 
situation  in  another  year.]  Our  current  solution  to  this  problem  is  the  Hessian-biased 
force  field 15  (HBFF)  which  combines  normal  mode  eigenstate  information  from  theory 
with  eigenvalue  information  from  experiment.  This  force  field  yields  accurate  vibrational 
frequencies,  and  we  have  been  systematically  determining  parameters  for  various  molecules. 

We  have  developed  a  new  program  (FFPARM)  that  is  used  in  conjunction  with  POLY- 
GRAF  to  determine  the  optimum  force  field  parameters  to  fit  a  set  of  experimental  data. 
The  input  data  consists  of  geometry,  crystal  cell  parameters,  elastic  constants,  vibrational 
frequencies,  phonon  frequencies,  or  the  full  Hessian  matrix  (from  quantum  mechanical  cal¬ 
culations).  These  data  can  come  from  theorj1-  or  experiment  (whichever  is  more  reliable) 
and  can  be  incomplete  (the  data  is  weighted  by  its  reliability).  FFPARM  then  finds  the 
optimum  force  field  parameters  (force  constants,  well  depths,  coupling  terms)  lor  valence 
and/or  vdw  terms  to  fit  the  input  data. 

A  number  of  projects  are  underway  (hydrocarbons,  fluorocarbons,  diamond,  silicon, 
germanium,  silicon-carbide,  amino  acids,  polyenes)  for  getting  accurate  force  fields  using 

FFPARM. 

Unfortunately  for  transition  metal  complexes,  there  is  generally  not  enough  experi¬ 
mental  data  for  this  approach  and  use  of  HF  wavefunction  is  not  adequate.  The  ability 
to  carry  out  PS-GVB-RCI  calculations  should  solve  this  problem  in  1992.  In  this  case  the 
calculated  constants  should  be  sufficiently  accurate  without  the  use  of  experimental  data. 

C.  van  der  Waals 

All  molecular  systems  lead  to  attractive  interactions  at  large  R  (dispersion  or  London 
forces)  and  repulsive  interactions  at  small  R  (Pauli  orthogonalization)  for  all  nonbonded 
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electrons.  Referred  to  as  van  der  Waals  interactions,  these  interactions  are  generally  as¬ 
sumed  to  be  pairwise  additive 

B.,h.  =  £  £?/"  (11) 

»>i 

and  describable  with  simple  exponential-6 

Ev{fw  =  Ae~BR:>  -  C/Rij  (12) 


or  Lennard-Jones  12-6  expressions. 

Our  long  range  goal  in  developing  force  fields  is  to  obtain  the  van  der  Waals  (disper¬ 
sion)  terms  directly  from  quantum  mechanical  calculations.  Unfortunately  the  weakness 
of  the  interactions,  the  complications  of  handling  the  dynamic  correlations  consistently 
and  completely,  and  the  need  for  very  complete  basis  sets  have  limited  reliable  calculations 
to  the  size  of  He2  and  Ne2.  MP2  or  EN2  level  wavefunctions  should  lead  to  acceptable 
accuracy  but  the  basis  sets  must  be  very  complete,  limiting  this  approach  to  quite  small 
molecules. 

As  a  result,  our  emphasis  is  on  determining  these  parameters  empirically  by  using 
new  simulation  procedures  to  predict  the  thermodynamic  and  mechanical  properties  of 
molecular  crystals  and  polymers  for  various  choices  of  the  van  der  Waals  parameters  and 
then  fitting  them  to  whatever  experimental  data  are  available.  We  have  now  developed 
programs  capable  of  calculating16,1' 

a.  the  full  set  (all  21)  of  elastic  constants  of  molecular  crystals  and  polymers. 

b.  the  vibrational  levels  and  phonon  dispersion  curves,  and 

c.  the  thermodynamic  parameters  (Cp,  S,  H,  G)  as  a  function  of  temperature. 

Although  the  experimental  data  are  incomplete,  we  find  that  b}'  calculating  all  such  prop¬ 
erties,  we  can  determine  the  optimal  van  der  Waals  parameters.17 

Thus,  to  determine  the  van  der  Waals  parameters  for  carbon  we  required  that  the 
lattice  constants  and  elastic  constants  be  accurately  described.18  Then,  to  obtain  the  van 
der  Waals  parameters  for  H,  we  required  that  the  crystal  structure,  heat  of  sublimation,  and 
lattice  frequencies  of  crystalline  polyethylene  be  accurately  described.19  This  was  tested 
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for  the  phonon  states  of  polyethylene  and  for  crystals  of  hydrocarbon  molecules  where  we 
calculated  accurate  lattice  parameters,  lattice  frequencies,  and  cohesive  energy.20  Similarly, 
new  van  der  Waals  parameters  have  been  obtained28  for  N,  0,  F,n  P,  S,  Cl,  Se,  Br,  Te, 
and  I,  and  others  are  in  progress. 

The  lowest  level  wavefunction  adequate  for  nonempirical  (ab  initio)  calculations  of 
van  der  Waals  parameters  is  second-order  perturbation  theory  (e.g.  MP2,  EN2).  With 
PS-MP2  and  PS-GVB-EN2  this  should  be  possible  for  systems  as  large  as  benzene  dimer. 
Such  studies  should  be  practical  by  mid  1992.  In  the  meantime,  we  are  using  FFPARM 
in  conjunction  with  POLYGRAF  (w'hich  calculates  phonon  bands,  elastic  constants,  and 
lattic  parameters)  to  determine  empirical  vdW  parameters  of  a  number  of  systems. 


D.  Metals  and  Alloys  (The  Interstitial  Electron  Model) 

Most  simulations  on  metallic  systems  have  used  two-body  force  fields,  and  the  results 
have  been  quite  poor.  Recent  advances  based  on  GVB  quantum  chemical  calculations  have 
led  to  the  Interstitial  Electron  Model  (IEM)  for  force  fields  of  metals  in  -which  the  electrons 
are  treated  as  classical  particles  on  the  same  footing  as  the  ionj!21  This  completely  new 
type  of  force  field  for  metals  will,  we  believe,  lead  to  accurate  simulations  for  metals. 

Based  on  generalized  valence  bond  calculations  of  small  metal  clusters  (up  to  14 
atoms),  we  learned  that  metallic  systems  lead  to  singly-occupied  orbitals  localized  at  in¬ 
terstitial  sites.  The  result  is  a  strong  counling  o:  electronic  states  with  nuclear  motion  so 


tnat  one  cannot  eescrioe  tne  sv 


stem  in  terms  of  ion  motionss  alone.  We  have  now  incor¬ 


porated  this  concept  of  intersiially  localized  orbitals  into  a  force  field  scheme  for  metals 
where  the  electrons  are  treated  as  classical  particles! 

This  force  field  is  then  written  as 


E  —  'y  '  (fini  +  'y  ^  <fie  +  y  ^  (fee ' 

ii 1  ic  ec' 

where  potentials  are  defined  f or  interactions  of  electrons  (ie  and  ee)  in  addition  to  the 
traditional  ii  interactions.  As  a  first  test  we  have  applied  this  scheme  to  the  metals  Ni, 
Cu,  Pd,  Ag,  Pt,  Au,  where  we  allow  twro  parameters  for  each  interaction  and  obtain  an 
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exact  fit  to  the  lattice  constant  and  the  three  elastic  constants,  and  an  excellent  fit  to  the 
phonon  dispersion  curves.  For  example,  for  Ni,  this  leads  to  an  exact  fit  to  the  experimental 
data  for  lattice  constants,  elastic  constants  and  phonon  dispersion  curves  (see  Figure  1). 
Note  in  particular  here  that  the  Cauchy  discrepancj'  —  C44  7^  0  is  correctly  described 
(not  possible  with  sums  of  ii  potentials  alone). 

We  are  now  expanding  these  IEM  calculations  to  predictions  of  surfaces  (for  Ni,  pre¬ 
liminary  calculations  indicate  that  the  surface  layer  relaxes  by  about  3%),  grain  boundaries, 
and  alloys.22 
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Figure  1. 
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With  GVB-PBC  we  will  be  able  to  get  the  IEM  parameters  directly  from  ab  initio 
calculations  (rather  than  experiment). 


E.  Force  Fields  for  Inorganic  and  Organometallic  Systems 

The  standard  force  fields  for  organic  and  biological  systems  (MM2,  MM3,  AMBER, 
CHARMM)  have  evolved  over  many  years  of  developement.  Yet  these  force  fields  deal 
primarily  with  just  four  atoms  (H,  C,  N,  0).  To  describe  inorganic  and  organometallic 
systems,  we  must  deal  with  the  other  100  atoms  of  the  periodic  table.  Our  strategy  to 
doing  this  is  to  use  generic  approaches,  where  the  number  of  independent  parameters  are 
few  and  vary  smoothly  as  a  function  of  rows  and  columns  in  the  periodic  table. 

Our  first  such  force  field12  is  DREIDING  (developed  in  collaboration  with  Dr.  S. 
Mayo  and  Dr.  B.  D.  Olafson  of  BioDesign)  which  handles  the  25  elements  in  the  B,  C, 

N,  0,  and  F  columns  and  the  Br,  C,  Si,  Ge,  and  Sn  rows  of  the  periodic  table.  Here 
there  are  two  parameters  per  atom  R°  and  8 where  R°  is  the  atomic  radius  and  0°  is 
the  equilibrium  angle  of  the  hydride  (e.g.  CH4,  NH3,  SH2,  etc.).  In  DREIDING,  the 
force  constants  are  determined  by  very  simple  rules  (e.g.,  K  —  TOO,  1400,  2100  kcal/A2  for 
single,  double,  and  triple  bonds)  that  are  independent  of  the  atom.  Despite  the  incredible 
simplicity,  this  force  field  gives  quite  accurate  geometries  for  a  wide  variety  of  structures. 
Thus,  using  a  random  selection  of  76  structures  from  the  Cambridge  database  (including 
such  groups  as  —NOT.  —  (S02)  — .  —  PO,  .  —  (POT'  — '  the  RMS  error  in  atom  oositions  is 

O. 225  A  with  average  errors  of  0.009  A  in  bonds.  0.57  degrees  in  bond  angles,  and  0.22 


degrees  in  torsions.12 

Professor  A.  K.  Rappe  (Colorado  State  University)  and  Goddard  have  collaborated  to 
develop  a  Generic  Force  Field  (GENFF)  suitable  for  inorganics  and  organometaiiics.23  In 
this  force  field  there  are  three  parameters  per  atom  rather  than  two!  The  new  parameter  is 
an  atomic  force  constant  Z°  defined  such  that  the  force  constant  for  bond  stretch  is  (13) 


where 


(13) 
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and  the  force  constants  for  a  bond  angle  i  —  j  —  l  is  given  by  a  similar  expression  involving 
Z°  and  Z°.  These  atomic  force  constants  Z°  vary  smoothly  along  rows  and  columns, 
allowing  them  to  be  interpolated  for  cases  without  data.  Thus  with  GENFF  one  can 
predict  structures  for  any  combination  of  elements  from  H  to  Lr  (element  103). 

Combined  with  QEq  which  can  predict  charges  for  any  element  from  H  through  Lr,  we 
have  a  good  start  for  simulating  any  inorganic  compound.  To  obtain  accurate  vibrational 
frequencies  will  no  doubt  require  additional  cross  terms  (angle-stretch,  etc.)  from  BHFF, 
however  GENFF  is  an  excellent  start  to  making  predictions  of  inorganic  structures. 


IV.  New  Developments  in  Molecular  Dynamics 
A.  Cell  Multipole  Method 

The  calculation  of  Coulomb  energies 


£ 


QiQi 

Rij 


(14) 


involves  N2/2  terms  for  N  particles.  Commercial  molecular  dynamics  software  all  involves 
processes  of  order  N2,  limiting  the  maximum  sizes  to  ~  5000  atoms.  For  many  problems 
it  is  necessary  to  do  an  atomic-level  simulation  on  systems  with  a  million  atoms.  Here 
the  ordinary  procedures  of  MD  are  just  not  practical.  We  have  recently  developed  the 
Cell  Multipole  Method  (CMM)  to  solve  this  problem.  The  f.rst  step  is  to  place  the  whole 
molecule  in  a  box  (200-500  A  on  a  side,  depending  upon  the  system'  which  is  divided 
into  eight  children  cels.  Each  child  cell  is  divided  into  cigh  grandchildren  cells,  etc.  until 
a  generation  is  reached  with  about  4  atoms  per  cell  (the  microcell).  For  a  million  atoms 
there  are  six  levels.  For  each  microcell  the  multipole  moments  (charge,  dipole,  quaarupole) 
are  calculated  and  these  are  shifted  and  added  to  obtain  the  mulitpole  momenets  of  the 
parent  cells,  grandparent  cells,  etc.  To  calculate  energies,  forces,  etc.  the  total  potential 
of  each  microcell  is  written  at 


Vtotal{R)  =  Vnear(R)  +  Vfar(R)  (15) 

where  the  interactions  with  all  atoms  in  the  26  neighboring  cells  are  included  in  Vneor 
but  the  interaction  with  the  1,000,000  -  100  more  distant  atoms  (Vyar)  are  all  evaluated 
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using  the  multipole  fields.  Here  the  farther  atoms  use  high  level  (parent)  cells  while  closer 
cehs  use  the  multipole  field  of  low  level  (child)  cells.  This  procedure  is  extremely  fast  and 
accurate  as  indicated  in  Figures  2  and  3. 
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Figure  i,.  Computational  costs  for  calculation  of  the  Coulomb  nonbond  inter¬ 
actions  for  a  series  of  poly(vinylidene  fluoride)  polymers  (N  =  1792, 
15360,  122880,  and  1013750  atoms). 
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B.  Canonial  Dynamics 

In  recent  years  a  number  of  methods  have  been  proposed  to  obtain  thermodynamic 
and  other  properties  directly  from  molecular  dynamics.  Nose,24  Parrinelloand  Rahman,2'" 
Nose  and  Klein,26  Andersen,27  etc.  have  proposed  modified  Lagrangians  such  that  the 
collection  of  momenta  and  coordinates  from  a  long  dynamics  calculation  satisfy  a  Boltz¬ 
mann  distribution  for  constant  pressure,  constant  temperature  conditions.  Vi  e  have  imple¬ 
mented  the  Nose  and  Rahman-Parrinello  techniques  for  periodic  systems  and  are  testing 
these  methods  for  various  properties,  including  thermal  expansion  coefficients  and  elastic 
constants.  These  methods  are  sensitive  to  internal  “mass”  parameters,  to  initialization, 
and  to  total  dynamics  time. 

C.  New  Developments  with  Monte  Carlo  Methods 

Our  protein/polymer  folding  project  has  led  to  programs  capable  of  combining 
Monte  Carlo  and  molecular  dynamics  techniques  for  rapid  sampling  of  configuration  space 
to  predict  optimum  conformations.  Currently  these  program  are  effective  at  folding  loops 
with  up  to  20  amino  acids.  This  is  useful  in  using  the  experimental  structure  for  one 
protein  to  predict  the  structure  of  homologous  proteins.  There  are  indications  that  Nose 
canonical  dynamics  will  lead  to  rapid  folding. 


V.  Molecular  Dynamics  Calculations  of  Friction 

In  carrying  out  calculations  of  bulk  or  surface  properties  of  solids,  v.v  use  periodic 
boundary  conditicns  to  extend  the  system  it.  three  dimensions  and  two  dimensions,  re¬ 
spectively.  In  this  case  the  calcuiational  effort  is  proportional  to  the  number  of  atoms  in 
the  unit  cell  and  no  artifical  surfaces  are  present.  Since  the  systems  are  infinite  in  three  or 
two  dimensions,  respec'iveiy.  one  must  be  very  careful  with  electrostatic  summations 


V" 


R.,\ 


and  dispersion  summations 


since  these  generally  converge  slowly.  The  general  approach  for  accelerating  the  conver¬ 
gence  of  such  sums  is  to  convert  each  sum  into  two  sums,  one  of  which  converges  quickly 
in  real  space  and  one  which  converges  quickly  in  reciprocal  space.  In  these  procedures 
there  is  a  parameter  77  which  controls  how  much  of  the  total  sum  goes  into  each  part.  We 
have  developed16  an  approach  for  dynamically  estimating  the  total  error  in  the  real  and 
reciprocal  space  parts  so  that  we  can  choose  the  7j  to  obtain  a  prespecified  level  of  accuracy 
in  the  minimum  calculation  time.  In  addition,  we  have  extended  the  formulae  to  include 
contributions  to  vibrational  frequencies  and  elastic  constants.16'25 

To  use  the  full  apparatus  of  PBC  for  problems  with  surfaces,  we  treat  the  system 
as  a  finite  slab  where  the  thickness  Z  is  large  enough  that  the  surfaces  do  not  interact,  see 
Figure  4. 


In  this  case  the  unit  cell  parameters  (of  the  two  dimensions  shown)  are  A  and  C.  The  space 
between  the  surfaces  can  be  a  vacuum,  a  gas  (H2,  02,  etc.),  a  liquid  (H20),  a  lubricant, 
or  a  second  solid  (see  Figure  5).  The  dynamics  must  be  carried  out  under  conditions  of 
constant  temperature  and  constant  pressure  in  order  to  simulate  a  system  in  a  normal 
constant  T,  p  environment.  In  addition,  we  can  apply  an  overall  stress,  say,  in  the  C 
direction,  and  calculate  the  lattice  parameters  and  elastic  properties  for  the  system  in 
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equilibrium  with  this  external  stress.  In  such  situations,  the  Lagrangian  L  =  T  —  V  is 
modified  to  include  additional  terms,  e.g., 

Lt  = - gkTlns  (1G) 

2  M,  *  K  1 

which  ensure  that  the  dynamics  trajectories  sample  phase  space  appropriately  for  a  canon¬ 
ical  distribution  for  a  system  in  equilibrium  with  a  heat  bath  at  temperature  T. 

To  stud}'  dissipation  and  internal  friction,  we  consider  an  oscillating  external  stress, 


Fx (0  =  Fx  cos u)t  (17) 

and  follow'  the  response  in  the  various  cell  coordinates 

=  H°af3  c°s(u;t  -  6a/3).  (18) 

The  phase  factor  tan  6ap  is  directly  related  to  the  energy  dissipation  and  can  be  mea¬ 
sured  experimentally  as  a  function  of  temperature  for  various  frequencies.  By  determining 
tz.n8ap  from  calculations  for  various  applied  loads  and  surface  treatments,  we  can  extract 
information  on  friction  for  various  surface  conditions.  Experimental  determinations  of 
tan  £  are  commonplace  for  polymers  where  they  provide  information  on  the  rigidity  of  the 
oolvmer  blend  as  a  function  of  temnerature.  Thev  have  no:  been  useful  for  stucvinc  sur- 
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oscillating  snear  (n),  -Jtnougn  tnese  calculations  will  interlace  only  indirectly  to  cur¬ 
rent  configurations  for  measuring  friction  it  should  be  possible  to  desicn  experiments  for 
studying  surfaces  such  as  mica.  There  may  be  some  gap  between  the  experimental  and 
theoretical  frequencies  since  our  current  programs  are  limited  by  practical  considerations 
(computer  time)  to  frequencies  above  10  GHz,  whereas  100  MHz  is  about  the  highest  used 
experimentally. 

For  diamond  and  silicon,  the  initial  focus  in  on  surfaces  passivated  with  H  or  F. 
Surfaces  where  all  surface  sites  are  capped  with  an  H  or  F  should  have  particularly  low- 


27 


friction.  Intuitively  the  fluorine  surface  might  be  expected  to  have  the  lowest  friction 
since  it  is  smooth  and  allows  little  in  the  way  of  lateral  interactions  due  to  surface-surface 
contact,  the  microscopic  origin  of  the  friction  in  such  systems  is  that  lateral  forces  on  the 
surface  atoms  due  to  the  contact  lead  to  bending  of  C-C-X  bonds  that  couple  with  the 
high  energy  branches  of  the  bulk  acoustic  phonon  modes.  This  coupling  should  be  stronger 
for  for  CCF  than  CCH  because  of  the  lower  frequencies,  and  hence  the  friction  may  be 
higher. 
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For  unpassivated  surfaces  having  broken  bonds,  additional-dissipation  effects  will 
quickly  increase  the  friction 


distorts  the  adjacent  bonds,  leading  to  coupling  into  various  pnonon  nodes  to  dissipate 
energy.  Our  calculations  should  be  able  to  quantify  these  effects. 


We  expect  that  the  distribution  of  phonon  modes  and  the  coupling  with  surface 
states  will  be  particularlj'  important  for  systems  such  as  diamond  and  silicon  and  have 
developed  general  procedures  for  calculating  these  values,  e.g.,  Figure  8  is  for  polyethylene, 
Figure  9  is  for  silicon. 

Metals 

We  see  very  different  mechanisms  for  friction  in  diamond  and  metallic  systems  and 
intend  to  carry  out  parallel  studies  on  clean  metal  systems.  Until  recently,  we  believed  that 
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it  was  impossible  to  carry  out  reliable  simulations  of  metal  systems  because  available  force 
fields  arc  quite  inadequate.  However,  the  interstitial  electron  model  described  in  Section 
III.D  proposed  a  genuine  breakthrough  here. 

We  are  now  expanding  these  calculations  to  predictions  of  surfaces  (for  Ni,  prelim¬ 
inary  calculations  indicates  that  the  surface  layer  relaxes  by  about  3%),  grain  boundaries, 
and  alloys.  We  believe  that  this  work  will  have  progressed  sufficiently  for  studies  of  friction 
in  about  12  months.  In  order  to  define  a  sj'stem  for  which  experiments  could  best  test  the 
theory,  we  plan  to  do  clean  Ag  and  Ni  passivated  with  S 
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Interaction  of  the  Ag  surfaces  with  H-terminated  diamond  should  be  a  good  case 
for  comparing  experiment  and  theory  since  Ag  could  not  react  with  the  H  of  diamond. 


Figure  8. 


29 


D 

O 


Figure  9.  Phonon  band  structure  of  silicon.  The  solid  line2  are  theoretical  predictions 
using  parameters  calculated  from  cluster  calculations.  The  triangles  and 
cricles  are  experimental  frequencies. 


VI.  Calculations  of  Friction 

In  order  to  gain  insight  into  the  extraction  from  microscopic  calculations  of  the 
quantities  (e.g.,  friction  coefficients)  useful  in  describing  macroscopic  phenomena,  we  are 
carrying  out  a  series  of  calculations  simulating  friction  between  two  diamond  surfaces.  The 
f.rst  calculations  examined  the  (ill)  surface  of  diamond  with  each  surface  atom  terminated 
by  a  C-H  bond.  In  order  to  eliminate  edge  effects,  we  treated  a  finite  thickness  slab  with 
periodic  boundary  conditions,  as  in  Figure  11. 
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The  lattice  constants  parallel  to  the  slab  were  restricted  to  remain  those  of  the  infinite 
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system;  however,  each  atom  was  allowed  to  move  an  equilibrium  position.  These  slabs 
were  then  stacked,  as  in  Figure  12a,  to  form  a  three-dimensional  infinite  system. 
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Figure  12. 
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The  spaces  between  the  slabs  and  the  relative  location  of  one  slab  with  respect  to  the  next 
were  optimized  to  obtain  the  overall  equilibrium  system.  This  system  represents  two  perfect 
diamond  (iii)  surfaces  in  equilibrium.  We  then  calculated  the  equilibrium  structures  as  each 
surface  is  displaced  with  respect  to  the  next,  while  simultaneously  imposing  an  external 
load  (10  GPa)  on  the  system.  This  was  carried  out  by  treating  the  lateral  displacement  as 
a  shear  in  the  fundamental  unit  ceil  of  the  three-dimensional  system,  as  indicated  in  Figure 
12b.  Thus,  for  each  value  of  shear  ( j3 ),  we  reoptimized  all  atomic  positions,  allowing  the 
surface  hydrogens  to  bend  out  of  the  wave  as  the  surfaces  were  displaced,  and  we  allowed 
the  distance  between  the  layers  to  reoptimize.  This  was  done  as  a  function  of  external 
stress.  The  results  are  presented  in  Figures  13  and  14.  Figure  13  shows  the  energy  surface 
(with  no  imposed  load)  superimposed  on  the  atom  position  of  the  lower  surface.  The 
minimum  position  has  the  C-H  bond  of  the  upper  surface  symmetrically  between  the  three 
C-H  bonds  of  the  lower  surface.  The  minimum  energy  path  for  displacement  has  the  C-H 
of  the  upper  surface  move  over  an  energy  barrier  of  0.14  kcal/mol  to  a  new  equilibrium 
position  with  the  C-H  of  the  upper  surface  above  a  second  layer  C  of  the  lower  system.  In 
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this  position,  the  energy  is  but  0.01  kcal  higher  than  the  original  equilibrium  position.  In 
the  minimum  energy  pathway,  the  CH  bond  of  the  upper  state  must  next  move  northeast 
toward  the  center  of  the  next  hexagon  and  then  east  to  the  next  second-layer  carbon  etc.  If 
the  surfaces  are  constrained  to  move  exactly  along  a  line,  the  lowest  energy  pathway  would 
be  one  that  crosses  over  the  sides  of  each  hexagon,  leading  to  a  barrier  of  0.47  kcal/mol. 
The  highest  energy  point  occurs  when  the  surfaces  are  displaced  so  that  the  lower  CH 
bond  points  at  the  CH  bond  of  the  upper  surface,  leading  to  an  energy  1.0  kcal  higher 
than  the  minimum  position. 
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Figure  14. 


Bat  Low  do  \v <:  gu  Ixoiu  this  optimum  energy  surface  to,  say,  a  coefficient  of  friction? 
.4s  one  approach,  consider  the  vertical  displacements  between  the  slabs  necesary  for  the 
system  to  move  along  the  minimum  energy  pathway  (see  Figure  14).  The  slabs  are  closest 
for  the  minimum  energy  position  where  the  CH  bond  of  the  top  slab  points  in  the  middle  of 
the  hexagon  of  the  lower  slab.  As  the  slabs  are  displaced  to  the  right  over  the  0.14  keal/mol 
barrier,  the  top  slab  must  move  away  by  0.061  A  and  then  it  comes  back  to  the  secondary 
minimum  (above  the  second  layer  C)  where  it  is  only  0.004  A  above  the  optimum  vertical 
displacement.  Thus,  as  the  system  is  displaced  along  the  minimum  energy  pathway,  the 
vertical  displacements  are  as  in  Figure  15. 

Vertical 

Displacement  0.061  0.061 

A 

0.  1.456  2.912 

Lateral  Displacement 

Figure  15. 


mere  is  an  external  loan  on  me  system,  we  must  do  work  arair.st  this  bond  to  disulace 
it  back  0.06  A  to  go  over  the  Lrst  hill,  but  the  bond  can  do  work  on  the  svstem  as  we  no 
corn  the  r.:..  again,  i:  tins  process  were  carrier,  out  with  :::■  enemy  transfer,  there  would 

and  is  men  available  to  cro  uo  tne  next  mil.  novever.  n  tne  cncrc v  is  mssmated  from, 
the  surface  as  it  goes  down  one  hill,  it  cannot  be  returned  to  climb  the  next  and  there  is 
triction.  To  get  an  idea  whether  there  is  time  for  such  dissipation,  consider  that  a  velocity 
Ox  2  meters/ sec  converts  to  0.02  A  per  picosec.  Thus,  at  this  velocity,  it  takes  35  picosec 
for  the  system  to  go  back  down  hill.  On  the  other  nand,  the  characteristic  time  for  energy 
transfer  from  the  C-C-H  bending  modes  at  the  surface  is  likely  to  be  ~1  psec.  Thus  the 
energy  should  be  dissipated. 

Given  the  description  in  Figure  15,  we  can  estimate  the  maximum  coefficient  of 
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friction  as  the  maximum  slope  of  the  displacement  curve  or 


f  =  7T 


0.061 

1.456 


=  0.132. 


It  is  interesting  that  this  value  is  in  the  ballpark  of  experimental  measurements  of  the 
coefficient  of  friction  in  diamond  by  Tabor  who  reports  f  =  0.1  for  diamond  (ill)-  This 
does  not  prove  our  simple  model.  The  experiments  on  diamond  were  carried  out  with 
a  metal  probe  in  a  normal  atmosphere.  We  now  know  that  a  well  polished  diamond 
surface  is  covered  with  C-H  bonds.  Hence  these  experiments  should  be  repeated  in  an 
inert  atmosphere  with  a  diamond  probe  so  that  one  does  not  have  to  worry  about  the 
surface  H  on  the  diamond  being  transferred  to  the  metal. 

It  is  also  not  clear  that  the  coefficient  of  friction  for  a  completely  fiat  surface  is  rel¬ 
evant  for  friction  studies  since  even  a  single  crystal  will  have  steps  and  other  imperfections 
that  may  keep  parts  of  the  surface  from  interacting  (Figure  16). 


1 
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That  is,  the  experiments  on  the  best  quality  single  surface  may  be  dominated  by  what 
happens  when  two  steps  come  together  rather  than  by  the  sliding  of  two  surfaces  with 
respect  to  each  other.  Of  course,  similar  theoretical  calculations  can  be  carried  out  at  the 
steps;  however,  there  are  new  questions.  Are  the  step  atoms  also  completel}'  saturated 
with  C-H  bonds?  We  could  assume  this  in  the  calculations,  but  this  should  be  tested  with 
appropriate  surface  science  experiments.  Using  these  theoretical  approaches,  one  should 
be  able  to  simulate  friction  and  wear  for  systems  such  as  Figure  16. 
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Over  the  next  year  we  will  be  carrying  out  similar  calculations  for  fluorine- 
terminated  surfaces  and  for  other  low  index  surfaces  of  diamond.  The  (110)  surface  should 
be  particularly  interesting  since  the  friction  might  be  very  anisotropic. 

VII.  Scanning  Tunneling  Microscopy 

A  new  tool  of  significant  potential  for  characterizing  the  surface  topography  of  tri¬ 
bologically  relevant  systems  is  scanning  tunneling  microscopy.  In  this  system,  a  small 
voltage,  say,  <f>  =  100  mV,  is  applied  between  a  probe  [prepared  so  as  to  obtain  extremely 
small  (atomic  dimensions)  points]  and  the  surface  and  the  tunneling  current  is  measured 
and  used  with  various  electronics  and  piezoelectric  elements  to  maintain  atomic-level  reso¬ 
lution.  The  probe  is  maintained  at  distances  of  5  to  10  A  and  the  tunneling  current  [which 
is  assumed  to  vary  exponentially  with  distance  (s)] 

IT(s)  =  exp  |  — A  sv^] 

is  measured.  In  one  mode  of  operation,  the  current  is  kept  constant  as  the  tip  is  displaced 
along  the  surface,  and  the  voltage  applied  to  the  piezoelectric  element  is  monitored  to 
obtain  the  effective  displacement  of  the  tip.  Assuming  that  the  work  function  of  the 
surface  6  is  independent  of  position,  one  then  gets  the  topography  of  the  surface.  The 
problem  here  is  that  at  the  atomic  level  the  effective  work  function  and  rate  of  electron 
transfer  does  depend  upon  the  nature  of  the  orbitals  describing  the  electrons. 

There  are,  however,  a  number  of  questions  about  tnese  experiments: 

1.  For  a  given  surface,  do  the  peaks  of  corrugation  occur  at  the  atomic  centers  or 
between  them? 

2.  For  heterogeneous  surfaces,  under  what  circumstances  can  STM  distinguish  atom 
types? 

3.  Under  what  operating  conditions  can  a  particular  adsorbate  be  imaged? 

4.  What  effects  does  tip  structure  have?  Can  we  tailor-make  more  effective  tips? 
An  STM  spectrum  was  taken  of  molybdenum  disulfide  by  the  Caltech  Group.1  A 

naive  interpretation  of  this  would  be  that  the  bright  spots  correspond  to  the  sulfurs,  since 
1  STM  data  courtesy  of  M.  Weimer  et  al.  of  Caltech  (1987). 
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these  are  farthest  from  the  surface  and  closest  to  the  tip.  The  intermediate  intensity  would 
be  the  molybdenums  which  are  one  layer  down,  and  the  lowest  intensities  would  be  the 
channels  that  would  go  entirely  through  the  system. 

The  way  we  are  currently  calculating  these  systems  is  to  represent  the  molybdenum 
disulfide  infinite  slab  using  clusters  as  shown  in  Figures  17  and  IS.  We  consider  that  this 
cluster  was  cut  from  the  infinite  system  and  that  wherever  there  is  a  broken  bond,  we  keep 
the  broken  bond  electrons  on  the  molybdenum  coupled  high-spin  (as  if  they  were  coupled 
to  the  down-spin  orbitals  of  the  missing  fragments  to  obtain  bond  pairs).  This  is  to  prevent 
these  electrons  from  rearranging  to  provide  stronger  bonding  within  the  complex  than  it 
would  have  in  the  infinite  system.  For  the  example  under  consideration  here,  there  are  30 
broken  bond  orbitals  and  28  electrons  that  should  be  distributed  among  these.  We  replace 
the  other  atoms  in  the  infinite  system  by  an  estimate  of  their  net  charges  (say,  molybdenum 
plus  two,  sulfur  minus  one)  and  then  calculate  this  finite  cluster  until  the  wavefunctions 
converge.  From  this  we  then  evaluate  the  net  charges  on  the  central  molybdenum  and 
sulfur  atoms  and  use  these  to  estimate  the  charges  on  the  atoms  that  are  missing  from 
the  system.  We  repeat  this  process  until  it  converges.  In  the  current  calculations,  the  tip 
is  modeled  with  a  Cs~  at  a  distance  R  from  the  surface  as  shown  in  Figure  19.  7  ~  each 
projection  of  cesium  above  the  surface  (Xh  )  and  for  each  distance  R  above  the  surface, 
we  calculate  two  wavefunctions.  One  denoted  TiTg  has  the  extra  electron  on  the  left  itt 


1';  )  and  the  other  wavefunction  denoted  'Ity  T'b  has  the  extra  electron  or.  th 


rate  or  electron  transier  is  given  dv 


o  — 

j  .  o 

Tlr  =  —  !Klr!‘  P- 
a 


where  ‘he  matrix  element  is  given  by 


Hlr  =  + 


and  H  is  the  full  many-electron  hamiltonian.  In  these  calculations  all  the  orbitals  of 
^a’I'b  overlap  all  the  orbitals  of  To  carry  out  the  calculation  for  Hlr,  we  use  the 


3S 


resonance-GVB  program  (developed  by  Voter  and  Goddard2)  that  allows  one  to  rapidly 
evalulate  such  matrix  elements  without  additional  approximations.  Results  from  these 
calculations  are  shown  in  Figures  20-22. 


2  Chem.  Phys.,  57,  253  (1981). 
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Figure  18. 
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Figure  19. 
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Figure  20. 
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Cs  located  above  Mo  —\ 


£  LECTIO  tJ  OU  cTs 

Cluster  HOMO  along  lines  through  probe  atom  (Cs) 


Figure  22.  Electron  on  Cs. 
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First,  in  Figure  20  we  consider  the  system  without  the  cesium;  that  is,  we  show 
just  the  density  of  the  orbital  on  the  MoS 2  that  gets  transferred  to  Cs+  but  without  the 
cps! um  present.  Tv,''co  densities  »re  non-monotor.ic  because  the  amplitude  of  the  orbital 
goes  to  zero  at  around  6  A  from  the  surface;  however,  overall,  these  are  higher  density 
above  the  sulfur  than  above  the  molybdenum. 

In  Figure  21  we  show  the  orbitals  for  the  'I'a'I'b  wavefunction,  that  is,  the  wave- 
function  with  the  electron  still  on  the  MoS 2.  Here  we  see  that  when  the  cesium  is  above 
the  molybdenum,  there  is  higher  density  from  about  5  A  out  than  when  the  cesium  is  above 
the  sulfur.  Figure  22  shows  the  case  where  the  electron  is  transferred  to  the  cesium,  that 
is  In  this  case,  the  amplitudes  are  essentially  the  same  from  about  4  A  out.  Thus 

the  electron  transfer  matrix  elements  Hlr  are  dominated  by  the  character  shown  in  Figure 
14. 

We  have  now  carried  out  calculations  of  the  relative  transfer  rates  for  these  two 
systems  and  find  that  indeed  when  the  tip  is  8.5  A  from  the  sulfur  (10  A  from  the  molyb¬ 
denum),  the  highest  electron  transfer  rate  is  when  the  tip  is  above  the  molybdenum !  That 
is,  the  measured  current  does  not  correspond  to  the  surface  topography.  We  are  in  the 
process  of  extending  these  results  in  several  ways. 

We  are  developing  an  approach  to  handle  an  infinite  slab  directly  so  that  one  does 
not  have  to  do  the  iterative  cluster  calculation  described  above  (which  complicates  inter¬ 
pretation).  In  addition,  we  are  using  the  current  approach  to  map  out  how  the  intensity 
changes  as  we  move  the  tip  around.  The  next  phase  of  this  work  will  be  to  carry  out 
similar  calculations  for  various  atoms  (hydrogen,  fluorine,  and  nothing)  sitting  on  various 
sites  above  diamond  and  silicon.  Because  the  STM  experiments  are  more  easily  carried 
out  on  silicon,  and  a  number  of  results  are  already  available,  we  will  examine  the  silicon 
sj’stem  first. 

Our  conclusion  from  these  studies  is  that  STM  has  a  great  deal  of  promise  for 
characterizing  tribological  systems;  however,  theory  will  be  critical  in  extracting  the  physics 
from  the  data. 
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17.  Simulation  of  Materials:  A-pplications  to  Biological,  Polymer,  and  Superconductor 
Systems,  Chemistry  Department  Seminar,  Yale  University,  New  Haven,  Connecti¬ 
cut,  28  February  1989.  (Presented  by  W.  Goddard.) 

18.  Simulation  of  Materials:  Applications  to  Biological,  Polymer,  and  Superconductor 
Systems,  General  Electric  Corporate  Research  and  Development  Center,  General 
Electric  Company,  Schenectady,  New  York,  1  March  1989.  (Presented  by  W.  God¬ 
dard.) 


19.  Simulation  of  Materials:  Applications  to  Biological,  Polymer,  and  Superconductor 
Systems,  1989  Joe  L.  Franklin  Memorial  Lecture,  Department  of  Chemistr}',  Rice 
University,  Houston,  Texas,  14  March  1989.  (Presented  by  W.  Goddard.) 

20.  Simulation  of  Biological  and  Materials  Systems,  Symposium  on  Computational 
Modeling  of  Molecular  Systems,  Divisions  of  Biological  Chemistry  and  Organic 
Chemistry  and  Subdivision  of  Theoreitcal  Chemistry,  197th  National  Meeting  of 
the  American  Chemical  Society,  Dallas,  Texas,  12  April  1989.  (Presented  by  W. 
Goddard.) 

21.  Simulations  on  Polymers  —  Mechanical,  Electrical  and  Structural  Properties,  Im¬ 
perial  Chemical  Industries  Molecular  Simulation  Symposium,  Wilton,  England,  20 
April  1989.  (Presented  by  W.  Goddard.) 

22.  Hessian-Biased  Force  Fields  for  Molecular  Mechanics,  West  Coast  Theoretical 
Chemistry  Conference,  IBM  Almaden  Research  Center,  San  Jose,  California,  11 
May  1989.  (Presented  by  S.  Dasgupta.) 

23.  Applications  of  Ab  initio  van  der  Waals  Potentials  of  H2—H2  and  CH4—H2  to  Large 
Molecular  Systems,  West  Coast  Theoretical  Chemistry  Conference,  IBM  Almaden 
Research  Center,  San  Jose,  California,  11  May  1989.  (Presented  by  J-G.  Lee.) 

24.  Structure  and  Bulk  Properties  of  Polymer  Crystals,  West  Coast  Theoretical  Chem¬ 
istry  Conference,  IBM  Almaden  Research  Center,  San  Jose.  California.  11  May 
1989.  (Presented  by  N.  Karasav.-a.) 


26. 


27. 


Molecular  Modeling  Applied  to  Problems  in  Polymers  and  Surface  Science.  Physical 
Chemistry  Department.  General  Motors  Research  Laboratories.  Warren,  Michigan. 
15  June  1989.  (Presented  by  W.  Goddard.) 

Atomic  Level  Simulation  of  Materials,  Workshop  in  Dedicated  High-Performance 
Systems  for  Visualization  and  Simulation,  Ramo  Auditorium,  Caltech,  Pasadena, 
California,  20  June  1989.  (Presented  by  W.  Goddard.) 

Simulation  of  Materials  with  Emphasis  on  Polymers,  Raychem,  Menlo  Park,  Cali¬ 
fornia,  28  June  1989.  (Presented  by  W.  Goddard.) 
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28.  Simulations  of  Materials:  An  Update,  General  Electric  Company,  Corporate 
Research  and  Development  Center,  Schnectady,  New  York,  August  9,  1989  . 
(presented  by  W.  Goddard) 

29.  Simulation  of  Polymers  and  Other  Materials,  Asahi  Glass  Reearch  Center, 
Yokohama,  Japan,  October  4, 1989  (presented  by  W.  Goddard) 

30.  Simulation  of  Polymers,  Superconductors,  and  Other  Materials,  Nippon  Steel, 
Kawasaki,  Japan  October  5, 1989  (presented  by  W.  Goddard) 

31.  Four  Lectures  on  Theory  and  Simulation  of  Chemical,  Biological,  and 
Materials  Systems,  Visiting  Professor  of  Chemistry,  Texas  A&M  University, 

College  Station,  Texas,  October  16-19, 1989  (presented  by  W.  Goddard) 

32.  Modeling  of  Chemical  and  Trihiological  Properties  of  Ceramic  System,  AFOSR 
Molecular  Dynamics  Contractors  Conference,  Session  on  Tribiological 
Surfaces,  Captiva  Island,  Florida,  November  1,  1989  (presented  by  W. 
Goddard) 

33.  Microscopic  Theoretical  Modeling  of  Ceramic  Interfaces,  Air  Force  Tribiology 

•4k 

Technical  Review,  Fairborn.  Ohio,  December  1,  1989  (presented  by  IV.  : 
Goddard)  r* 

34.  Simulation  of  Materials:  Applications  to  Polymers ,  Metals,  and  Catalysts,  BP 
America  Research,  Cleveland,  Ohio,  March  12,  1990  (presented  by  W. 
Goddard) 

35.  Tutorial  on  Molecular  Modeling,  Exxon  Chemical  Company  Chemicals 
Research  Meeting,  Buck  Hill,  Pennsylvania,  May  1,  1990  (presented  by  W. 
Goddard) 

36.  Simulations  of  Surfaces  and  Interfaces,  at  Alcoa  Laboratories  Technical 
Symposium  on  Chemistry  and  Physics  of  Surface  and  Interfaces,  Alcoa  Center, 
Pennsylvania,  August  20, 1990  (presented  by  W.  Goddard) 
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